home *** CD-ROM | disk | FTP | other *** search
/ Disc to the Future 2 / Disc to the Future Part II Programmer's Reference (Wayzata Technology)(6013)(1992).bin / MAC / THINKC / 4_0 / ORBMECHD / KEPLER.C < prev    next >
Text File  |  1991-02-18  |  8KB  |  358 lines
  1. /*****************************************************************
  2.     Kepler's equation solution algorithms from Battin, Mathematics
  3.     and Methods of Astrodynamics, pp192-6.
  4.     
  5.     8 August 1990
  6. *****************************************************************/
  7.  
  8. #include     "orbmech.h"
  9. #include    "SANE.h"
  10. #include    "math.h"
  11. #define        KEPLERREC_SIZE        50
  12.  
  13. /* globals */
  14.  
  15. extern    DialogPtr        KeplerDia;
  16. extern    EventRecord        gEvent;
  17. extern    TEHandle        TEH;
  18. extern    decform            gdecform;
  19. extern    int                dirty;                    
  20.  
  21. struct    ParamRec
  22.     {
  23.     extended    value;
  24.     Str255        name;
  25.     };
  26.     
  27. extern     struct    ParamRec    gParameter;
  28.  
  29. struct    KeplerRec
  30.     {
  31.     extended         f_1, f_2, time;
  32.     extended         a, e;  
  33.     } ;
  34.     
  35. typedef        struct    KeplerRec    *KeplerPtr;
  36.     
  37.  
  38. /*****************************************/
  39.  
  40. Kepler()
  41. {
  42.     Handle        inputH, itemHandle;
  43.     Str255        f1Str, f2Str, ftStr, aStr, eStr;
  44.     int            itemType;
  45.     Rect        itemRect;
  46.     Boolean        dont_panic = 1, miFlag = 0, auFlag = 0, degFlag = 0;
  47.     decimal        result;
  48.     
  49.     ControlHandle    kmSecBut, miHrBut, radBut, timeBut;
  50.  
  51.     struct        KeplerRec    *kdata;
  52.     
  53.     Ptr            temp;
  54.  
  55.     Str255        q, w, e, r, t, y, u, i;    
  56.     
  57.     temp = NewPtr( KEPLERREC_SIZE );
  58.     
  59.     kdata = ( KeplerPtr) temp;
  60.  
  61.     GetDItem( KeplerDia, F_TIME, &itemType, &inputH, &itemRect);
  62.     GetIText( inputH, &ftStr );
  63.     Pstr_to_extended( &ftStr, &kdata->time, &dont_panic );
  64.     
  65.     GetDItem( KeplerDia, F_1, &itemType, &inputH, &itemRect);
  66.     GetIText( inputH, &f1Str );
  67.     Pstr_to_extended( &f1Str, &kdata->f_1, &dont_panic );
  68.     
  69.     GetDItem( KeplerDia, F_2, &itemType, &inputH, &itemRect);
  70.     GetIText( inputH, &f2Str );
  71.     Pstr_to_extended( &f2Str, &kdata->f_2, &dont_panic );
  72.     
  73.     GetDItem( KeplerDia, A_IN, &itemType, &inputH, &itemRect);
  74.     GetIText( inputH, &aStr );
  75.     Pstr_to_extended( &aStr, &kdata->a, &dont_panic );
  76.     
  77.     GetDItem( KeplerDia, E_IN, &itemType, &inputH, &itemRect);
  78.     GetIText( inputH, &eStr );
  79.     Pstr_to_extended( &eStr, &kdata->e, &dont_panic );
  80.     
  81.     GetDItem ( KeplerDia, KM_SEC_O, &itemType, &kmSecBut, &itemRect);
  82.     GetDItem ( KeplerDia, MI_HR_O, &itemType, &miHrBut, &itemRect);
  83.     GetDItem ( KeplerDia, RAD_BUT, &itemType, &radBut, &itemRect);
  84.     GetDItem ( KeplerDia, TIME_UNK, &itemType, &timeBut, &itemRect);
  85.     
  86.     if ( dont_panic )
  87.     {
  88.     
  89.     if ( ! GetCtlValue ( kmSecBut ) )
  90.         if ( GetCtlValue ( miHrBut ) ) {
  91.             kdata->a *=  MI_CONV;
  92.             kdata->time *=  HR_CONV;
  93.             miFlag = 1;
  94.         }
  95.         else {
  96.             kdata->a *=  AU_CONV;
  97.             kdata->time *=  YR_CONV;
  98.             auFlag = 1;
  99.         }    
  100.     
  101.     if ( ! GetCtlValue ( radBut ) )  {
  102.         kdata->f_1 /=  DEG_CONV;
  103.         kdata->f_2 /=  DEG_CONV;
  104.         degFlag = 1;
  105.     }
  106.     
  107.     if ( GetCtlValue ( timeBut ) )  {   /*  start time solution */
  108.         
  109.         pStrCopy("\p\rTime between two points ", q);
  110.         pStrCopy("\p\r     start point:    ", w);
  111.         pStrCopy("\p\r     end point:      ", e);
  112.         pStrCopy("\p\r     a:              ", r);
  113.         pStrCopy("\p\r     e:              ", t);
  114.         pStrCopy("\p\r     primary:        ", y);
  115.         pStrCopy("\p\r\rTime:                ", u);
  116.         pStrCopy("\p\r\r", i);
  117.         
  118.         TEInsert( q+1, *q, TEH );
  119.         TEInsert( w+1, *w, TEH );
  120.         TEInsert( f1Str+1, *f1Str, TEH );
  121.         TEInsert( e+1, *e, TEH );
  122.         TEInsert( f2Str+1, *f2Str, TEH );
  123.         TEInsert( r+1, *r, TEH );
  124.         TEInsert( aStr+1, *aStr, TEH );
  125.         TEInsert( t+1, *t, TEH );
  126.         TEInsert( eStr+1, *eStr, TEH );
  127.         TEInsert( y+1, *y, TEH );
  128.         TEInsert( gParameter.name+1, *gParameter.name, TEH );
  129.         
  130.         KeplerFindTime( kdata );
  131.         
  132.         if ( miFlag ) 
  133.             kdata->time /=  HR_CONV;
  134.         
  135.         if ( auFlag )  
  136.             kdata->time /=  YR_CONV;
  137.         
  138.         if ( kdata->time > 9999.0  || kdata->time <  .000001 ) 
  139.             gdecform.style = FLOATDECIMAL;
  140.         else
  141.             gdecform.style = FIXEDDECIMAL;
  142.         
  143.         num2dec( &gdecform, kdata->time, &result );    
  144.         dec2str( &gdecform, &result, ftStr );
  145.         
  146.         TEInsert( u+1, *u, TEH );
  147.         TEInsert( ftStr+1, *ftStr, TEH );
  148.         TEInsert( i+1, *i, TEH );
  149.  
  150.         dirty = 1;
  151.         
  152.         GetDItem( KeplerDia, F_TIME, &itemType, &itemHandle, &itemRect);
  153.         SetIText( itemHandle, ftStr );
  154.  
  155.         ShowSelect();
  156.     }
  157.     else  {         /**** find position block ***/
  158.         
  159.         pStrCopy("\p\rPosition after elapsed time", q);
  160.         pStrCopy("\p\r     start point:    ", w);
  161.         pStrCopy("\p\r     time:           ", e);
  162.         pStrCopy("\p\r     a:              ", r);
  163.         pStrCopy("\p\r     e:              ", t);
  164.         pStrCopy("\p\r     primary:        ", y);
  165.         pStrCopy("\p\r\rFinal position:      ", u);
  166.         pStrCopy("\p\r\r", i);
  167.         
  168.         TEInsert( q+1, *q, TEH );
  169.         TEInsert( w+1, *w, TEH );
  170.         TEInsert( f1Str+1, *f1Str, TEH );
  171.         TEInsert( e+1, *e, TEH );
  172.         TEInsert( ftStr+1, *ftStr, TEH );
  173.         TEInsert( r+1, *r, TEH );
  174.         TEInsert( aStr+1, *aStr, TEH );
  175.         TEInsert( t+1, *t, TEH );
  176.         TEInsert( eStr+1, *eStr, TEH );
  177.         TEInsert( y+1, *y, TEH );
  178.         TEInsert( gParameter.name+1, *gParameter.name, TEH );
  179.         
  180.         KeplerFindPosition( kdata );
  181.         
  182.         if ( degFlag ) 
  183.             kdata->f_2 *=  DEG_CONV;
  184.  
  185.         if ( kdata->f_2 > 9999.0  || kdata->time <  .000001 ) 
  186.             gdecform.style = FLOATDECIMAL;
  187.         else
  188.             gdecform.style = FIXEDDECIMAL;
  189.         
  190.         num2dec( &gdecform, kdata->f_2, &result );    
  191.         dec2str( &gdecform, &result, f2Str );
  192.         
  193.         TEInsert( u+1, *u, TEH );
  194.         TEInsert( f2Str+1, *f2Str, TEH );
  195.         TEInsert( i+1, *i, TEH );
  196.  
  197.         dirty = 1;
  198.         
  199.         GetDItem( KeplerDia, F_2, &itemType, &itemHandle, &itemRect);
  200.         SetIText( itemHandle, f2Str );
  201.  
  202.         ShowSelect();
  203.     }
  204.         
  205.     }
  206.  
  207. }
  208.  
  209. /*****************************/
  210.  
  211. extended    trueToEccentric( f, e )
  212. extended    *f, *e;
  213. {
  214.     extended    ecc;
  215.     
  216.     ecc = 2 * atan( sqrt( (1 - *e) / (1 + *e) ) * tan( 0.5 * *f ) );
  217.     
  218.     return (ecc);
  219. }
  220.  
  221. /*****************************/
  222.  
  223. extended    eccentricToTrue( ecc, e )
  224. extended    *ecc, *e;
  225. {
  226.     extended    f;
  227.     
  228.     f = 2 * atan( sqrt( (1 + *e) / (1 - *e) ) * tan( 0.5 * *ecc ) );
  229.     
  230.     return (f);
  231. }
  232.  
  233. /*****************************/
  234.  
  235. KeplerFindTime( kdata )
  236. struct    KeplerRec    *kdata;
  237. {
  238.     extended    f1, f2, k, a, period, ecc1, ecc2, m1, m2;
  239.     extended    t1, t2, e, temptime, tf1, tf2;
  240.     
  241.     f1 = kdata->f_1;
  242.     f2 = kdata->f_2;
  243.     tf1 = kdata->f_1;
  244.     tf2 = kdata->f_2;
  245.     a = kdata->a;
  246.     e = kdata->e;
  247.     
  248.     k = sqrt( a * a * a / gParameter.value );
  249.     
  250.     period = 2 * PI_VALUE * k;
  251.  
  252.     if ( tf1 > ( PI_VALUE - 0.00001 ) && tf1 < ( PI_VALUE + 0.00001 ) )  
  253.         t1 = period / 2.0;
  254.         
  255.     else  {
  256.         
  257.         if ( tf1 >= ( PI_VALUE - 0.00001 )  )
  258.             tf1 = 2 * PI_VALUE - tf1;
  259.         
  260.         ecc1 = trueToEccentric( &tf1, &e );
  261.         
  262.         t1 = k * ( ecc1 - e * sin( ecc1 ) );
  263.     }
  264.  
  265.     if ( tf2 > ( PI_VALUE - 0.00001 ) && tf2 < ( PI_VALUE + 0.00001 ) )  
  266.         t2 = period / 2.0;
  267.         
  268.     else  {
  269.         
  270.         if ( tf2 >= ( PI_VALUE - 0.00001 )  )
  271.             tf2 = 2 * PI_VALUE - tf2;
  272.         
  273.         ecc2 = trueToEccentric( &tf2, &e );
  274.         
  275.         t2 = k * ( ecc2 - e * sin( ecc2 ) );
  276.     }
  277.     
  278.     temptime = 0;
  279.     
  280.     if ( ( f1 > f2 ) || ( f2 > PI_VALUE && f1 < PI_VALUE ) )
  281.         temptime += period;
  282.     
  283.     if ( f1 < PI_VALUE )
  284.         t1 *= -1.0;
  285.     
  286.     if ( f2 > PI_VALUE )
  287.         t2 *= -1.0;
  288.     
  289.     temptime += t1 + t2;
  290.     
  291.     kdata->time = temptime;
  292.     
  293. }
  294.  
  295. /***************************************/
  296.  
  297. KeplerFindPosition( kdata )
  298. struct     KeplerRec    *kdata;
  299. {
  300.     extended    k, a, period, ecc1, m1, t1, ttot, time, f1;
  301.     extended    f2, ecc2, m2, t2, old_ecc, e;
  302.     
  303.     int            converged;
  304.     
  305.     f1 = kdata->f_1;
  306.     
  307.     time = kdata->time;
  308.     
  309.     a = kdata->a;
  310.     
  311.     e = kdata->e;
  312.     
  313.     k = sqrt( a * a * a / gParameter.value );
  314.     
  315.     period = 2 * PI_VALUE * k;
  316.  
  317.     if ( f1 > ( PI_VALUE - 0.00001 ) && f1 < ( PI_VALUE + 0.00001 ) )  
  318.         t1 = period / 2.0;
  319.         
  320.     else  {
  321.         
  322.         if ( f1 >= ( PI_VALUE - 0.00001 )  )
  323.             f1 = 2 * PI_VALUE - f1;
  324.         
  325.         ecc1 = trueToEccentric( &f1, &e );
  326.         
  327.         t1 = k * ( ecc1 - e * sin( ecc1 ) );
  328.     }
  329.     
  330.     ttot = t1 + time;
  331.     
  332.     while ( ttot > period )
  333.         ttot -= period;
  334.         
  335.     converged = 0;
  336.     
  337.     old_ecc = 0;
  338.     
  339.     ecc2 = 0;
  340.     
  341.     m2 = ttot / k;
  342.     
  343.     while ( ! converged )  {
  344.     
  345.         ecc2 = m2 + e * sin( ecc2 );
  346.         
  347.         converged = checkConvergence( ecc2, old_ecc );
  348.         
  349.         old_ecc = ecc2;
  350.     }
  351.     
  352.     f2 = eccentricToTrue( &ecc2, &e );
  353.     
  354.     if ( f2 < 0.0 )
  355.         f2 += 2 * PI_VALUE;
  356.     
  357.     kdata->f_2 = f2;
  358. }